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Abstract 

A unitarized BFKL equation incorporating shadowing and antishadowing cor- 
rections of the gluon recombination is proposed. This equation reduces to the 
Balitsky-Kovchegov evolution equation near the saturation limit. We find that the 
antishadowing effects have a sizeable influence on the gluon distribution function in 
the preasymptotic regime. 
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1 Introduction 

The growth of cross sections with gluon splitting according to the DGLAP [1] and 
BFKL [2] evolution equations would violate the unitarity. Therefore, the corrections of the 
higher order QCD, which shadow the growth of parton densities and lead to an eventual 
saturation of parton densities, become a focus of intensive study in recent years. In this 
respect, the GLR-MQ (by Gribov, Levin Ryskin in [3] and Mueller and Qiu in [4]) and BK 
(by Balitsky and Kovchegov in [5]) equations are broadly regarded as two solvable models 
to compute the shadowing corrections to the DGLAP and BFKL equations, respectively. 

One of the key problems in restoring unitarity is the origin of the negative corrections. 
In the viewpoint of elementary QCD interaction, the suppression of the gluon splitting 
comes from its inverse process-the gluon recombination. The negative screening effects 
in the recombination processes originally occur in the interferant cut-diagrams of the 
recombination amplitudes [3,6]. For computing the contributions from the interference 
processes, the AGK cutting rules [7] were used in the derivation of the GLR-MQ equation. 

On the other hand, one of us (Zhu) disputed the above mentioned applications of the 
AGK cutting rules in the GLR-MQ equation, and used the TOPT-cutting rules based 
on the time ordered perturbation theory (TOPT) instead of the AGK cutting rules to 
expose the relations among various cut diagrams in a recombination process [6]. Thus, we 
can completely compute the contributions of the gluon recombination processes. Conse- 
quently, a modified DGLAP equation was proposed [6,8]. A remarkable property of this 
equation is that the positive antishadowing and negative shadowing components in the 
nonlinear evolution equation are naturally separated. As a result, the corrections from 
gluon recombination at small x depend not only on the size of gluon density at this value 
of x, but also on the shape of gluon density in the region [x/2,x]. Thus, the shadowing 
effects in the evolution process will be weakened by the antishadowing effects if the gluon 



distribution has a steeper form. 

The antishadowing effects always coexist with the shadowing effects in the QCD recom- 
bination processes as a general conclusion of the momentum conservation [9]. Therefore, 
similar antishadowing effects should exist in any unitarized BFKL equation. In this work 
we try to find the antishadowing effects in the unitarized BFKL equation. The total 
momentum carried by gluons is not changed as a result of the gluon fusion but simply 
redistributed to different x regions. The depleted momentum due to the fusion is heaped 
up step by step toward a larger x direction, which leads to an enhanced density of higher 
momentum gluons (i.e. the antishadowing effects). We shall show that the antishadowing 
effects have a sizeable influence on the gluon distribution function in the preasymptotic 
regime. 

The paper is organized as follows. In Section 2 we derive a unitarized BFKL equation 
with the gluon recombination, which contains the contributions of the shadowing and 
antishadowing effects. The numerical analysis of our equation are presented in Section 3. 
In Section 4 we show that our equation reduces to the BK equation near the saturation 
limit. We also compare our equation with a modified BK equation in this Section. 



2 The evolution equation incorporating shadowing 

and antishadowing effects 

We develop a unified framework to construct the evolution equations for both the 
integrated and unintegrated gluon distribution functions. We begin from a deep inelastic 
scattering process, where the unintegrated gluon distribution is measured. In the kr- 
factorization scheme, the cross section is decomposed into 



da(probe*P -> k'X) 
fk 2 x \ 

= f(xi,kl) <8> /C =o, — , a s ® da(probe*k — > fc') 
\«i ^1 / 

= Af(x,k 2 )0da(probe*k^k'), (1) 

which contains the perturbative evolution kernel /C, the nonperturbative unintegrated 
gluon distribution function / and the pro6e*-parton cross section da(probe*k — > k'). For 
simplicity, we take the fixed QCD coupling in this work. According to the scale-invariant 
parton picture of the renormalization group [10], we regard Af(x, k 2 ) as the increment of 
the distribution f(xi,k%) when it evolves from (xi,k%) to (x, k 2 ). Thus, the connection 
between the two functions f(xi,kj) and fix,]? 2 ) via Eq.(l) is 



f(x,k 2 ) = f(x u kj) + Af(x,k 2 ) 

=f ^r^/cf^f,^/^,^, ( 2 ) 

4? mi „ h Jx x 1 v^i x i / 

The relation of the unintegrated gluon distribution with the integrated gluon distri- 
bution is 

rQ 2 dk 2 rQ 2 dk 2 

G(x,Q 2 ) = xg(x,Q 2 )=J 2 ^xf(x,k 2 ) = 2 ^F(x,k 2 ), (3) 

—rain — —rain — 



or 



ft 7 n 2 d 9{x,Q 2 ) 

f(x,k ) = Q dQ2 Q2= fc2 . (4) 



In the evolution along the transverse momentum (or along the longitudinal momen- 
tum), we differentiate Eq.(2) with respect to A; 2 (or to x) and get 

df(x,k 2 ) 
dk 2 

2 



+ : % / 7 v -^r J f(^kD, (5) 
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or 



df(x,k 2 ) 
dx 

* dk 2 1 /A; 2 x \ 2 
^ 72'— > a « J\ x ^h) 



k_ ^ X\ \^k_\ ^i 



respectively. Generally, the resummation solution is hard to obtain from these two equa- 
tions. However, at the LL(k 2 )A (or at the LL{l/x)A) the evolution kernel /C in Eq.(5) (or 
in Eq.(6)) is only the function of the longitudinal variable (or of the transverse variables) 
and the second term on the right-hand side of Eq.(5) (or of Eq.(6)) vanishes. In this case, 
the resummation becomes possible. For example, using Eq.(3) we write 

A^,Q 2 )= f $-Af(x,k 2 ) 

f Q 2 dk 2 ft 2 dkl ridx! fx \ 

= /„ 7f , "72"/ I^dglap — ,a s }J[Xi,k l ) 

4Ln k 2 Jk\ min k 2 Jx Xl \ Xl J 

= [, TT / —^dglap (—,®s) g(xi,k 2 ), (7) 
4L„ k Jx xi \xi J 



and 



g(x,Q 2 )=g(x 1 ,k 2 ) + Ag(x,Q 2 ). (8) 
On the other hand, at DLL(Q 2 )A we have 

dx 1 a s N c dxx 

i^DGLAP = • (yj 

X\ 7T X 

Thus, from Eqs.(5) and (7) one can get the DGLAP equation at small x for the gluon 
distribution 

2 dg(x,Q 2 ) f l dx 1 fx \ 2 

= a I N 1 ^d_x 1 x_ 1 

7T Jx X\ X 

Now let us consider the corrections of the gluon recombination to the DGLAP equa- 
tion. The one step evolution containing the gluon recombination is illustrated in Fig. 1(a), 
where the initial gluons from the two legs are resummed using the DGLAP equation. 
Four different kinds of the evolution kernels for JCmd-dglap are listed in Fig.2. Similar 
to Eqs.(7) and (8) we derive the modified DGLAP equation from 

G(x,Q 2 ) = G(x l ,Q 2 1 )+AG(x,Q 2 ) 
= G(x 1 ,Q?)-2 f f 1/: 1 ^-IC M ^ DGLAP (-,a s )G^(x u Ql) 

J Qlmin Ql Jx X l X l \X! J 

+ /T ^ [ 1/2 —-£M D - D GLAp(-,a s )GW(x 1 ,Q 2 1 ), (11) 
jQlmiu Qi Jx / 2 x i x i \xi J 

where a power suppression factor 1/Q\ has been extracted from the evolution kernel. The 
positive and negative terms are from the contributions of the diagrams Fig.2 (a) and (b), 
respectively. The TOPT calculations [8] give 



j^MD-DGLAP d x l 
Xi 

a 2 N 2 (2 Xl - x)(72x\ - A8x\x + lA<dx\x 2 - llQx^ + 29x 4 ) 

~ y w^t xjx~ 1 

— 18a s ^3T— • ( 12 ) 

The gluon correlation function G^ 2 ^ is a generalization of the gluon distribution beyond 
the leading twist. It is usually modelled as the product of two gluon distributions [6]. For 
example, 

G^ 2 \x,Q 2 )=R G G 2 (x,Q 2 ), (13) 

where R G = l/(irR 2 ) is a correlation coefficient with the dimension [L~ 2 ], R is the effective 
correlation length of two recombination gluons and we take R = SGeV' 1 (the radius of a 
nucleon). Setting Eqs.(12) and (13) to Eq.(ll), we obtain the modified DGLAP equation 
combining DGLAP dynamics at small x 

dG{x,Q 2 ) 
dlnQ 2 

a s N c f 1 dx 1 2 36a 2 N 2 f l/2 dx, 2 2 
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nQ 2 R 2 N 2 -llx,2~ G {XuQ) - (14) 



It is interesting to compare this small- a; version of the modified DGLAP equation with 
the GLR-MQ equation, which is [4] 



dG(x,Q 2 ) 



dlnQ 2 



_a s N c [ 1 dx 1 2 36a 2 N 2 C dx 1 2 2 

where 

G (2) (^g 2 ) = ^G 2 (x,g 2 ), (16) 

is assumed. Although both of the two equations have been hoped to describe the correc- 
tions of the gluon recombination to the linear DGLAP equation at the DLLA, they have 
different forms due to the following reasons: 

(1) The GLR-MQ equation is derived basing on the following two works. One is 
the idea that the negative shadowing effects arise from the interference processes of the 
gluon recombination. Gribov, Levin and Ryskin used the AGK cutting rules to show 
that the contributions from the real cut diagrams (see Fig. 2 (a)) and the interferant cut 
diagrams (see Fig.2(b)) only differ in the numerical weights, which is 2 and -4, respectively. 
Thus, the net effects of the gluon recombination can be simply calculated by multiplying 
the real diagrams by a negative weight. Obviously, the antishadowing effects in this 
approach are completely cancelled by the shadowing effects and the resulting evolution 
equation violates the momentum conservation. Another work is the calculation of the 
recombination functions at the DLLA in a covariant perturbation theory by Mueller and 
Qiu [4], where a special treatment was used to remove the infrared (IR) singularities in the 
gluonic twist-4 coefficient functions. The gluon recombination functions in the GLR-MQ 
equation are separated out from the above mentioned IR-safe coefficient functions. 

(2) In the derivation of the modified DGLAP equation, the TOPT was first used 
to establish the relations among all the relevant cut diagrams in Fig. 2 [6]. We showed 
that the shadowing and antishadowing effects share the same recombination function but 
occupy different kinematic regions. On the other hand, the contributions of the two 
virtual diagrams Fig. 2(c) and (d) cancel against each other. Thus, the net effects depend 



not only on the local value of the gluon distribution at the observed point, but also on 
the shape of the gluon distribution when the Bjorken variable goes from x to x/2. In 
consequence, the shadowing effects in the evolution process will be obviously weakened 
by the antishadowing effects if the distribution is steeper [11]. We found that the same 
TOPT framework allows us to calculate the recombination functions including quarks 
and gluons in the whole region of x at the LL(Q 2 )A. Consequently, the momentum 
conservation is restored in the modified DGLAP equation [6]. We also showed in [8] that 
the gluon recombination functions can be reasonably separated from the twist-4 coefficient 
functions in the TOPT approach without the above mentioned special treatment, though 
the results of the two methods only differ by a constant at the DLL A. 

In next step, we consider the recombination of two unintegrated gluon distribution 
functions, which obey the BFKL equation. The cut diagram corresponding to Fig. 1(a) is 
Fig. 1(b), where the dashed box implies a new recombination kernel of two BFKL solutions. 
The initial gluons couple with the gluonic dipole and it means that two initial legs evolve 
according to the BFKL dynamics. One can expect that JCmd-bfkl is more complicated 
than JCmd-dglap- As an approximate model, we use kernel JCmd-dglap to replace 
the kernel ICmd-bfkl m Fig. 1(b). In the concrete, the contributions of two correlated 
unintegrated distribution functions to the measured (integrated) distribution G via 
the recombination processes are 

AG(x,Q 2 ) = -2 f % r ———-JCmd—dglap (-,<*.) F^\x u k 2 ) 

J QLn t Jx Xl Xl \Xl J 

fQ 2 dk 2 r 1 / 2 dxi x „ fx \ , 21 2 

+ / -fr \ ——ICMD-DGLAp[—,a s )F {2 >(x 1 ,k 2 ). 17 

jQ 2 m in k J x/2 Xi Xi VXi J 

Thus, we have 



AF(,,a = Q 2SAG(X ' Q2; 



ag 2 

2 r 1 / 2 dxi x ( x \ (2 \, ,2s 

j2 iCmd-dglap — ,a s j F K >(xi,k ) 

K_ Jx X\ X\ \X\ 

1 f 1 / 2 dxi x .„ fx 



1 r 1 / 2 dx\ x fx \ / 9 x. , 9 , , . 

+72/ — — /CjifD-cGLAPI— ,a a F^(a;i,fc 2 ). 18 

X 1 X 1 \Xi J 

Using Eqs.(6) and (12) we obtain the corrections to the evolution of the unintegrated 
gluon distribution along small x direction 



dF( X ,k 2 ) 2 N 2 J_ F (2)/ fc 2x , lg 2 ,1^(2)^ ^ / 1Q v 

Combining with the BFKL equation, we obtain a unitarized BFKL equation 



dF(x,k 2 ) 
dx 

a s N c e r°° dk! 2 \ F(x,tf 2 ) ~ F(x,k 2 ) F{x,ti 



7T 4:r . /-- /2 I l/. y2 - - /• J: h . /& 4 + 4^' 4 



-n ? ?. T7 



nk 2 R 2 N 2 — 1 v ' - y nk 2 R 2 N 2 - 1 
where similar to Eq. (13) we define 



F^ 2 \x,k 2 ) = -^F 2 (x,k 2 )- (21) 

Eq.(20) is our unitarized BFKL equation for the unintegrated gluon distribution. Since 
the gluon distribution becomes flatter near the saturation limit, in this range we can take 
the approximation 



in Eq.(20), and get 



F^,k 2 ^F(x,k 2 ), 



(22) 



dF(x,k 2 ) 
dx 

a s N c k 2 dk! 2 { F(x,g 2 ) - F(x,k 2 ) F(x,k 2 ) 
4^ ~F \ |i' 2 -i 2 | + 
18a? N 2 



F 2 {x,k 2 ). (23) 



nk 2 R 2 - 1 



3 Numerical analysis 

We make some numerical calculations to illustrate the antishadowing effects. For 
simplicity, we fix the coupling constant to be a s = 0.3. Eq.(20) is the corrections of the 
gluon recombination to the BFKL evolution equation. Therefore, we use a parameter 
form of the BFKL solution as the input distribution at x = 1CT 2 [12] 

where Xbfkl = 12o; s /(7rln2) and A" ~ 32.1o; s . The parameter /c 2 = lGeV 2 is of nonper- 
turbative origin and the normalization constant (3 ~ 0.01. On the other hand, we take 
££ m = 0.01Ge\/ 2 in Eq.(20). 

We assume that Eq.(20) begins to work from the BFKL-region. In this case, the 
distribution F(xo, fc 2 ) has a steeper x-dependence at the starting point x of the evolution 
where the shadowing and antishadowing effects almost cancel out. Thus, we compute the 
value of F(xi/2,k 2 ) at the i-th step of the evolution from F(xi,k 2 ) using the BFKL 
equation, i.e., 

F( Xl ,k 2 ) BFK ^ path F (|, A; 2 ) = Fbfkl (f ,£ 2 ) . (25) 

The numerical results show a faster divergence of the distribution F(x,^) due to the an- 
tishadowing effects that cancel out or even outweigh the shadowing effects in the typical 
BFKL-solution. Of course, we can not expect such infinite growth of the gluons at small 
x. This unphysical result appears because the BFKL dynamics are not constrained by the 
energy-momentum conservation, which must be added from the outside of the evolution 
equation. According to the experiences in the solutions of the modified DGLAP equa- 
tion, the net antishadowing effects in a small x region always accompany the stronger net 



shadowing effects in the smaller x region [11]. One can expect that the evolution asymp- 
totically approaches the BK dynamics under the action of the net shadowing effects. 
Therefore, we modify the program (25) to 

p ( x * h A - F ( x% lcA 4. Fbfkl (iiz!) ~ Fshadowm9 (fjj) 

V 2~' ) ~ hshadowm 3 ) + irj F -ri F + l ' ( ' 

where Fshadowing(xi/2, k 2 ) indicates that the evolution from %i to Xi/2 is controlled by 
Eq.(23). The parameter rf implies the different velocities approaching the BK dynamics. 
There is a limit value of rf lin ~ 0.001 and the solution of Eq.(20) becomes divergent when 
rf < rf nin . As an example, we take rf = rf nin and divide the evolution into the 100-steps 
on |Alogx| = I. 

In Fig. 3(a) we plot the unintegrated gluon distribution function (solid curves) as a 
function of x for two different values of k 2 and where r]^ in = 0.001. The possible solutions 
of Eq.(20) should lie between the solid and point curves since rf > rf^. For compar- 
ison, the results of two calculations based on the BFKL equation (dashed curves) and 
Eq.(23) without the antishadowing effects (point curves) are listed. One can find that 
Eq.(20) keeps the BFKL-behavior in a larger evolution region than Eq.(23) due to the 
antishadowing effects. 

The /c 2 -dependence of the unintegrated gluon distribution function with rf^ = 0.001 
is given in Fig.3(b). The saturation of distribution is usually defined as a limit form, 
which is insensitive to a; or fc 2 . We have not observed such saturation phenomenon in a 
broad kinematical range, although the obvious suppressions of the net shadowing effects 
can be found in Fig. 3. The reason is that the factor 1/k 2 in both the Eqs.(20) and (23) 
suppresses the contributions of the nonlinear terms at larger values of A; 2 . 

We give the integrated gluon density using Eq.(3) in Fig. 4. We still have not observed 
the saturation limit on the Q 2 -dependence of the distribution G(x, Q 2 ) in Fig. 4(b). It 



is interesting that this conclusion is qualitatively consistent with our previous work [11], 
where the corrections of the same gluon recombination kernel to the DGLAP evolution 
equation are considered. Of course, to obtain more practical predictions about the gluon 
distribution at small x, further corrections should be considered, for example, the cor- 
rections from the running a s , the BFKL dynamics in the IR-region h/ 2 -C lGeV 2 , and 
in particular, the choice of the input distribution. However, the important distinction 
between the antishadowing effects and shadowing effects as demonstrated in Fig.4 should 
remain after those corrections are incorporated. 



4 Discussions 



We try to compare the nonlinear evolution equation (20) with the BK equation, which 
is originally written in the transverse coordinate space for the scattering amplitude. There 
are different definitions of the scattering amplitude in the QCD evolution equation. In 
this work, we note that the BFKL equation for the scattering amplitude N(k?, x) has the 
following form [13], 



dN(k 2 ,x) 
-x- 



dx 

a s N c p dk! 2 \ k! 2 N{k! 2 , x) - k 2 N(k 2 , x) k 2 N{k 2 ,x) \ 

one can find that F(x,k 2 ) and N(k 2 ,x) differ by a simple scale translation F(x,k 2 ) ~ 
k?N(k 2 ,x). Thus, we define the scattering amplitude 

N(k 2 ,x) = ^F(x,k 2 ), (28) 

where we choose such constants on the right-hand side of equation that the following 
resulting equation is consistent with the BK equation at the saturation limit. 
Using Eq.(28) we rewrite Eq.(20) as 



dN(k 2 ,x) 
-x- 



dx 

, 1 , A',. p dk! 2 I k! 2 N{k! 2 , x) - k 2 N{k 2 , x) k 2 N{k 2 , x 



n Jk'^ k' 2 I \k' 2 -k 2 



yj k 4 + 4k! 4 

_2^N 2 (k 2 ,x) + ^N 2 (k 2 , ^ . (29) 

7T 7T V 2/ 

Interestingly, this equation is consistent with the BK equation (in the impact parameter- 
independent case) in momentum space near the saturation limit due to 



(30) 



For the cylindrically symmetric solution, we use the following transformation according 
to [14] 



/d 2 v 
— — exp(—ik ■ r) 
2,71 

r°° dr 



N(r, x) 



= / — J (kr)N(r,x), (31) 
Jo r 

where Jo is the Bessel function, one can find that Eq.(29) has the following form in 
transverse coordinate space 



-x- 



d N(r m x) 
dx 



OsNc 

2tt 2 



d 2 r 



-60 



-bcLcO 



[N(r bc , x) + N(r c0 , x) - N(r b0 , x) 



-2N(r bc , x)N(r c0 , x) + N (r bc) |) N [r c0 , |)], 



(32) 



which reduces to the BK equation (in the impact parameter-independent case) at the 
saturation limit 



— x 



dN(r b0 ,x) 
dx 



a s N c 
2tt 2 



, [AT( rfec) x ) + N(r c0 , x) - N(r b0 , x) - N(r bc , x)N(r c0 , x)] . 

L bc 31c0 



(33) 



We take the Golec-Biernat-W-usthoff model [15] as the input amplitude, i.e., 



N(r, x ) = 1 — exp 



r 2 Q' 2 



(34) 



where Q' s ~ lAGeV relates to the saturation scale. The computing program of Eq.(32) 



is 



N[k,-)=N BK (k,-) + 171 n_ 7]N + 1 > (35) 



where N BK (k 2 ,Xi/2) indicates that the evolution from Xj to Xj/2 is according to Eq.(33). 
The minimum parameter rj^ in ~ 0.01, which is larger than rj^in since the nonlinear effects 
in Eq.(32) are stronger than that in Eq.(20). 

Fig. 5 shows the solutions (solid curves) of Eq.(32) with r]^ in ~ 0.01. For comparison, 
we give the solutions of the BK equation (point curves) and the linear part of Eq.(32) 
(dashed curves). The possible solutions of Eq.(32) should lie between the solid and point 
curves since r/ N > rj^ in . Fig.6 is such a possible solution with r/ N = 0.1. One can find that 
the antishadowing effects lead to quite different behaviors between the amplitude N(r, x) 
and the distribution F(x, A; 2 ). The antishadowing effects on the scattering amplitude 
sharpen the transition form of the amplitude N(r, x) and even violate the unitarity due 
to iV > 1 if we take the definition Eq.(28). The reasons for the stronger antishadowing 
effects are that (1) we use a normalized definition about the amplitude in Eq.(28) so that 
the coefficients of the nonlinear terms are the same magnitude as that of the linear terms 
in Eq.(32); (2) the suppression factor 1/fc 2 is absorbed by the scattering amplitude due 
to Eq.(28). 

It is interesting that the momentum compensation for the shadowing effects was 
also discussed using a modified nonlinear evolution equation in Ref. [16]. In the impact 
parameter-independent case this equation reads 



_ x dN(rbo,x) _ a s N c f r 2 



dx 2tt 2 J c licdo 

^2N(r bc ,x) - N(r bc ,x)N(r c0 ,x) - ^[2N(r bc , x) - N(r bc , x)N(r c0 , x)]j , (36) 

where Y = ln(l/x) and the last term is the leading order DGLAP corrections. For 
comparison, we rewrite our Eq.(32) as 



— X 



dN(r m ,x) 



dx 



OsNc 

2tt 2 



/ 




{N(r bc , x) + N(r c0 , x) - N(r b0 , x) - 2N(r bc , x)N(r c0 , x) 



OsNc 

2tt 2 



+N(r bc ,x)N(r c0 ,x) - - — [N(r bc , x)N(r c0 , x)} j 
d 2 Lc^\- {N(r bc , x) + N(r c0 , x) - N(r b0 , x) - N(r bc , x)N(r c0 , x) 



L bc LcQ 



(37) 



Although Eqs.(36) and (37) have similar nonlinear antishadowing terms (only differ by 
a factor 1/2), one can find that the main difference between these two equations is that 
the same momentum compensation mechanism acts on both linear and nonlinear terms 
in Eq.(36), while the corresponding corrections in Eq.(37) are from the virtual and re- 
combination processes, respectively. 

Although we have discussed the contributions of the antishadowing effects to the 
unitarized BFKL equation, we note, however, that the following problems remain to be 
solved: (i) As we have emphasized, the replacement of the recombination functions for two 
BFKL-solutions with Kmd-dglap in this work is an approximate method. The unitarized 
BFKL equation near the saturation range should take a new form ; (ii) How to determine 
the parameter rf in Eq.(26)? Perhaps, the most reliable test of the antishadowing effects 
would be a comparison to the nuclear antishadowing effects, which has been observed in 
the EMC effect [17]. However, the kinematic region of the nuclear antishadowing effects is 
about 0.05 < x < 0.3, where we should combine the DGLAP dynamics in Eq.(20) and it is 
beyond the scope of the present paper; (iii) The QCD evolution should be performed with 
running a s in the numerical procedure. These problems will be our following subjects. 

In conclusion, we presented the corrections of the gluon recombination to the BFKL 
equation and they lead to a new unitarized nonlinear evolution equation, which incorpo- 



rates both shadowing and antishadowing effects. The new equation reduces to the BK 
equation near the saturation limit. The numerical solutions of the equation show that the 
antishadowing effects have a sizeable influence on the gluon distribution function in the 
preasymptotic regime. 
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Figure Captions 

Fig. 1 The cut diagram containing 4 2 gluon recombination kernel (the dashed 
box) in one step QCD evolution: (a) for the modified DGLAP equation and (b) for the 
modified BFKL equation, where two legs before fusion are resummed using the DGLAP 
and BFKL equations in (a) and (b), respectively. The dark area implies the probe. 

Fig. 2 The cut diagrams of the gluon recombination kernels for (a) the real processes 
that yield the antishadowing effects; (b) the interference processes that yield the shad- 
owing effects in the modified DGLAP equation, (c) and (d) are the corresponding virtual 
diagrams. 

Fig. 3 (a) The unintegrated gluon distribution function F(x,k?) as the function of x 
for different values of k 2 . The solid-, point- and dashed-curves are the solutions of Eq.(20) 
with ?7^j n = 0.001, Eq.(23) and BFKL equation, respectively; (b) Similar to (a) but as 
the function of k 2 for different values of x. The possible solutions of Eq.(20) should lie 
between the solid and point curves. 

Fig. 4 Similar to Fig.3 but for the gluon distribution function G(x, Q 2 ). 

Fig. 5 (a) The normalized scattering amplitude N(r, x) as the function of x for different 
values of r; (b) Similar to (a) but as the function of r for different values of x. The solid-, 
point- and dashed-curves are the solutions of Eq.(32) with 7]^ in = 0.01, the BK equation 
(33) and the linear part of Eq.(32), respectively. The possible solutions of Eq.(32) should 
lie between the solid and point curves. 

Fig. 6 A possible solution of Eq.(32) with r/ N = 0.1 (a) as the function of x, and (b) 
as the function of r. 
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